% This function runs the second step of the SR approach with specific
% options for what is a allowed to switch per equation
function [theta2,AVarTheta2Step2] = step2SR_threshold_macro_spec(...
    AVarTheta1Step1,outputStep1,bandwidthT,econAct,...
    switch_h0,switch_hx,numBootStep2)

dispOn         = 1;
allowDeSaling  = 1;           % allow for descaling of the var_u to ensure omegeHat is positive definite
thresholdLevel = 0;
[theta2,deScalingX,deScalingZ,hx_1,hx_2] = Var1Threshold_GMM_closedForm_macro_spec(...
    AVarTheta1Step1.VarFactorsRobust,...
    AVarTheta1Step1.autoCov10Factors,...
    AVarTheta1Step1.autoCov01Factors,...
    outputStep1.factors,outputStep1.macros',...
    econAct,thresholdLevel,allowDeSaling,switch_h0,switch_hx);
if dispOn == 1
    if deScalingX ~= 1
        disp(['The measurement error was descalted by ', num2str(deScalingX), ' in the estimated P dynamics for x'])
    end
    if deScalingZ ~= 1
        disp(['The measurement error was descalted by ', num2str(deScalingZ), ' in the estimated P dynamics for z'])
    end
end

% % Getting starting values
% if any(switch_h0 == 1)
%     h0RegimeOn = 1;
% else
%     h0RegimeOn = 0;
% end
% if any(any(switch_hx == 1))
%     hxRegimeOn = 1;
% else
%     hxRegimeOn = 0;
% end
% allowDeSaling  = 1;           % allow for descaling of the var_u to ensure omegeHat is positive definite
% thresholdLevel = 0;
% [theta2] = Var1Threshold_GMM_closedForm_macro(...
%     AVarTheta1Step1.VarFactorsRobust,...
%     AVarTheta1Step1.autoCov10Factors,...
%     AVarTheta1Step1.autoCov01Factors,...
%     outputStep1.factors,outputStep1.macros',...
%     econAct,thresholdLevel,allowDeSaling,h0RegimeOn,hxRegimeOn);
% % Removing elements in theta2Start to match switch_h0 and switch_hx
% nx = length(switch_h0);
% for i=1:nx
%     if switch_h0(i,1) == 0
%         theta2 = rmfield(theta2,['h0',num2str(i),'_2']);
%     end
% end
% for i=1:nx
%     for j=1:nx
%         if switch_hx(i,j) == 0
%             theta2 = rmfield(theta2,['hx',num2str(i),num2str(j),'_2']);
%         end
%     end
% end
% selectTheta2 = fieldnames(theta2);
% theta2Values = struct2array(theta2)';
% 
% % Optimization
% if all(switch_h0 ==1) && all(all(switch_hx))
%     theta2ValuesOpt = theta2Values;
% else
%     options         = optimset('Display','iter','MaxIter',5D4,'MaxFunEvals',5D4,'TolFun',1D-10,'TolX',1D-8);
%     theta2ValuesOpt = lsqnonlin(@objectFuncStep2Threshold_macro_spec,theta2Values,...
%         [],[],options,selectTheta2,AVarTheta1Step1.VarFactorsRobust,...
%         AVarTheta1Step1.autoCov10Factors,AVarTheta1Step1.autoCov01Factors,...
%         outputStep1.factors,outputStep1.macros',econAct,thresholdLevel,switch_h0,switch_hx);
% end
% theta2 = values2struct(theta2ValuesOpt,selectTheta2);
% % Adding common elements across regimes
% for i=1:nx
%     if switch_h0(i,1) == 0
%         name1 = ['h0',num2str(i),'_1'];
%         name2 = ['h0',num2str(i),'_2'];
%         theta2.(name2) = theta2.(name1);
%     end
% end
% for i=1:nx
%     for j=1:nx
%         if switch_hx(i,j) == 0
%             name1 = ['hx',num2str(i),num2str(j),'_1'];
%             name2 = ['hx',num2str(i),num2str(j),'_2'];
%             theta2.(name2) = theta2.(name1);
%         end
%     end
% end

% Ensure that the threshold model is stationary
[delta,theta2] = runInduceStatVAR_threshold_macro_spec(theta2,...
    AVarTheta1Step1.VarFactorsRobust,...
    AVarTheta1Step1.autoCov10Factors,...
    AVarTheta1Step1.autoCov01Factors,...
    outputStep1.factors,outputStep1.macros',econAct,thresholdLevel,switch_h0,switch_hx);

if numBootStep2 > 0
    seedNum = 1;
    numBootStep1And3 = outputStep1.numBootStep1And3;
    disp('RUNNING BIAS-ADJUSTMENT OF P-DYNAMICS')
    outBiasCorrect  = biasCorrectTheta2_threshold_macro_spec(theta2,AVarTheta1Step1.VarFactorsRobust,...
         AVarTheta1Step1.autoCov10Factors,...
        AVarTheta1Step1.autoCov01Factors,...
        outputStep1.factors,outputStep1.macros',econAct,thresholdLevel,...
        bandwidthT,switch_h0,switch_hx,allowDeSaling,numBootStep2,seedNum,numBootStep1And3);
    theta2 = outBiasCorrect.theta2Adj;
    delta  = outBiasCorrect.delta;
end


% The standard errors for theta2 using asymptotical inference
AVarTheta2Step2 = getAVarTheta2_threshold_macro_spec(theta2,AVarTheta1Step1.VarFactorsRobust,...
    AVarTheta1Step1.autoCov10Factors,AVarTheta1Step1.autoCov01Factors,...
    outputStep1.factors,outputStep1.macros',econAct,bandwidthT,thresholdLevel,switch_h0,switch_hx);
AVarTheta2Step2.namesFull    = fieldnames(theta2);

% Additional output
AVarTheta2Step2.delta = delta;
AVarTheta2Step2.thresholdLevel = thresholdLevel;
AVarTheta2Step2.switch_h0      = switch_h0;
AVarTheta2Step2.switch_hx      = switch_hx;
end

